#!/usr/bin/env python3
# coding: utf-8
"""
This module contains classes that represents the internal coordinates of a
molecular systems.
All classes inherit from the `InternalCoord` class.
"""
from pymatgen import Element
__all__ = ["Bond", "Angle", "Dihedral", "Improper"]
class InternalCoord:
""" A general class for internal coordinates """
def __init__(self, atoms, value, frc_cste, elements):
"""
Store internal coordinates data:
Args:
atoms (tuple): list of atom indexes(int(iat), int(jat))
value (float): equilibrium value of the coordinates
frc_cste (float): force constant value
elements (tuple): list of elements, either string or pymatgen.Element
"""
self.atoms = tuple(atoms)
self.value = float(value)
self.frc_cste = float(frc_cste)
self.elements = tuple(elements)
def as_dict(self):
"""
Return a serializable dict to export the internal coordinate to a json
file. Look at monty.json and MSNoable class for the future. The approach
is the one followed by pymatgen.
"""
return {"atoms": self.atoms, "value": self.value,
"frc_cste": self.frc_cste,
"elements": [el.specie.symbol for el in self.elements]}
@classmethod
def from_dict(cls, d):
"""
Return an InternalCoord object from a dictionnary.
Args:
d (dict): A dict representation of the internal coordinate such as
{"atoms": (1, 2), "value": 1.32, "frc_cste": 324.7, "elements": ["C", "N"]}
Returns:
InternalCoord object
"""
defaults = set(("atoms", "value", "frc_cste", "elements"))
args = set(d.keys())
if args != defaults:
# missing arguments
missing = defaults - args
if missing:
raise KeyError("Data are missing:\n%s" % "\n".join(missing))
# additionnal arguments
additionnal = args - defaults
if additionnal:
print("WARNING: The following data will not be considered\n%s" %
"\n".join(additionnal))
d["elements"] = [Element(el) for el in d["elements"]]
return cls(**{k: d[k] for k in defaults})
def __repr__(self):
line = self.__class__.__name__
line += "("
for k, v in vars(self).items():
line += "%s=%s, " % (k, str(v))
return line[:-2] + ")"
def __eq__(self, other):
"""
assume that two internal coordinates that are defined from the same atom
indexes are the same (equal)
"""
return self.atoms == other.atoms
def __hash__(self):
""" hash method based on the hash of the tuple of atom indexes """
return hash(tuple(self.atoms))
[docs]class Bond(InternalCoord):
""" A class that represents a bond """
def __init__(self, *args, **kwargs):
"""
Store internal coordinates data:
Args:
atoms (list): list of atom indexes
value (float): equilibrium value of the coordinates
frc_cste (float): force constant value
elements (list): list of elements, either string or pymatgen.Element
"""
super().__init__(*args, **kwargs)
[docs]class Angle(InternalCoord):
""" A class that represents an angle """
def __init__(self, *args, **kwargs):
"""
Store internal coordinates data:
Args:
atoms (list): list of atom indexes
value (float): equilibrium value of the coordinates
frc_cste (float): force constant value
elements (list): list of elements, either string or pymatgen.Element
"""
super().__init__(*args, **kwargs)
[docs]class Dihedral(InternalCoord):
""" A class that represents an dihedral angle """
known_pot = ("harmonic", "parabolic")
def __init__(self, atoms, value, frc_cste, elements, periodicity=None,
potential="parabolic"):
"""
Store internal coordinates data:
Args:
atoms (list): list of atom indexes
value (float): equilibrium value of the coordinates
frc_cste (float): force constant value
elements (list): list of elements, either string or pymatgen.Element
periodicity (int): periodicity in case of harmonic or sinusoidal potential
potential (str): functional form of the potential energy
"""
super().__init__(atoms=atoms, value=value, frc_cste=frc_cste, elements=elements)
self.periodicity = periodicity
if potential not in self.known_pot:
print("WARNING: potential is %s, and should be " % potential, self.known_pot)
print(" Be sure you know what you are doing.")
self.potential = potential
[docs] def as_dict(self):
"""
Return a serializable dict to export the internal coordinate to a json
file. Look at monty.json and MSNoable class for the future. The approach
is the one followed by pymatgen.
"""
d = super().as_dict()
d["periodicity"] = self.periodicity
d["potential"] = self.potential
return d
[docs] @classmethod
def from_dict(cls, d):
"""
Return a dihedral object from a dictionnary.
Args:
d (dict): A dict representation of the dihedral angle such as
{"atoms": (1, 2, 3, 4), "value": 120.32, "frc_cste": 4.7,
"elements": ["H", "C", "N", "H"], "potential": "harmonic",
"periodicity": 3}
Returns:
Dihedral object
"""
args = set(d.keys())
defaults = set(("atoms", "value", "frc_cste", "elements"))
if args != defaults:
# missing arguments
missing = defaults - args
if missing:
raise KeyError("Data are missing:\n%s" % "\n".join(missing))
# additionnal arguments
additionnal = (defaults | {"periodicity", "potential"}) - defaults
if additionnal:
print("WARNING: The following data will not be considered\n%s" %
"\n".join(additionnal))
d["elements"] = [Element(el) for el in d["elements"]]
if "periodicity" in d:
defaults.add("periodicity")
if "potential" in d:
defaults.add("potential")
return cls(**{k: d[k] for k in defaults})
[docs] def compute_periodicity(self, nnj, nnk):
"""
Compute the periodicity of a dihedral angle according to the procedure
proposed by Nilson et al. (Acta Cryst 2003 D59, 274 - 289).
The dihedral angle is supposed to be a torsion around atom j and atom k
which are bounded together and correspond to the second and the third
atoms of the self.atoms field of the dihedral.
::
i l
\ /
j -- k
Args:
nvoisj (int): number of neighbors of atom j
nvoisk (int): number of neighbors of atom k
Returns:
The periodicity (int)
"""
if self.periodicity is not None:
print("WARNING: the current periodicity value will be overwrite." +
" The current value is ", self.periodicity)
if self.potential != "harmonic":
print("WARNING: You ask for the calculation of the periodicity but" +
"the potential is not harmonic.")
if nnj > 4 or nnk > 4:
print(self, nnj, nnk)
raise ValueError("Maybe there is a problem with your dihedrals" +
" (a metal atom?). If this is not the case," +
" contact the authors of Mammoth, please.")
elif nnj == nnk:
if nnj == 2:
p = 1
elif nnj == 3:
p = 2
elif nnj == 4:
p = 3
else:
if max(nnj, nnk) == 2:
p = 1
else:
p = 6
self.periodicity = p
[docs]class Improper(InternalCoord):
""" A class that represents an improper dihedral angle """
def __init__(self, *args, **kwargs):
"""
Store internal coordinates data:
Args:
atoms (list): list of atom indexes
value (float): equilibrium value of the coordinates
frc_cste (float): force constant value
elements (list): list of elements, either string or pymatgen.Element
"""
super().__init__(*args, **kwargs)